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ABSTRACT 

Planetesimals embedded in a protoplanetary disc are stirred by gravitational torques exerted 
by density fluctuations in the surrounding turbulence. In particular, planetesimals in a disc 
supporting fully developed magneto-rotational turbulence are readily excited to velocity dis- 
persions above the threshold for catastrophic disruption, halting planet formation. We aim to 
examine the stirring of planetesimals lying instead in a magnetically-decoupled midplane dead 
zone, stirred only by spiral density waves propagating out of the disc's magnetically-coupled 
turbulent surface layers. We extend previous studies to include a wider range of disc models, 
and explore the effects of varying the disc column density and external magnetic field strength. 
We measure the stochastic torques on swarms of test particles in 3D resistive-MHD stratified 
shearing-box calculations with ionisation by stellar X-rays, cosmic rays, and recombination on 
dust grains. The strength of the stirring is found to be independent of the gas surface density, 
which is contrary to the increase with disc mass expected from a simple linear wave picture. 
The discrepancy arises from the shearing out of density waves as they propagate into the dead 
zone, resulting in density structures near the midplane that exert weaker stochastic torques on 
average. We provide a simple analytic fit to our numerically obtained torque amplitudes that 
accounts for this effect. The stirring on the other hand depends sensitively on the net vertical 
magnetic flux, up to a saturation level above which magnetic forces dominate in the turbulent 
layers. For the majority of our models, the equilibrium planetesimal velocity dispersions lie 
between the thresholds for disrupting strong and weak aggregates, suggesting that collision 
outcomes will depend on material properties. However, discs with relatively weak magnetic 
fields yield reduced stirring, and their dead zones provide safe-havens even for the weakest 
planetesimals against collisional destruction. 
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1 INTRODUCTION 

The formation of rocky earth-like planets, and the rock-ice cores 
of gas giant planets, are believed to involve assembly from a pop- 
ulation of planetesimals that may range in size from 1 to 100 km. 
The existence of planets, asteroids and Kuiper-belt objects in our 
own solar system covering a wide-range in semimajor axes, in ad- 
dition to the broad diversity in orbital properties of extrasolar plan- 
ets (e.g. [Mayor & Queloz|1995[|Marois et aI.|20T0{ |Lissauer et al.| 
|201 1| > and debris discs around main-sequence stars ( Wyatt 2008[l, 
indicates that planetesimals form at orbital distances ranging from 
a few tenths of an au out to more than 100 au from the central star 
While there remain many open questions about how planetesimals 
form within a protoplanetary disc, a number of plausible models 
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and scenarios have been suggested that exhibit differing degrees of 
sensitivity to the physical conditions within the disc. 

Classical models of planetesimal formation (see e.g., |Wei-| 
[denschilling & Cuzzi|1993] > involve an incremental process which 
starts from micron-sized dust grains. These grains grow via mutual 
collisions and sticking and settle toward the midplane of the proto- 
planetary disc (PPD). Assuming a laminar nebula a few times heav- 
ier than the minimum required to form the solar system ( |Hayashi| 
1 1 98 1| >, [Weidensch illing (200()i predicts a growth time scale equiv- 
alent to a few xlO'' yr at 5 au. The main obstacle to the incremental 
formation scenario is the rapid inward migration of growing boul- 
ders once they reach sizes around one metre. Moreover, as demon- 
strated by [Brauer et al.| |2008|, differential radial migration of a 
population of planetesimals of varying sizes can lead to high rela- 
tive velocities implying collisional fragmentation. 

In an attempt to avoid these difficulties, alternative forma- 
tion scenarios have been proposed, which try to jump this metre- 
sized barrier by involving mechanisms of rapid planetesimal for- 
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mation. These alternative scenarios go back to the classic idea of 
|Goldreich & Ward| ( |1973 '), who proposed formation of planetes- 
imals through gravitational instability in a dense layer of solids 
resulting from vertical settling of dust grains. Alternative scenar- 
ios that invoke rapid planetesimal formation include the concentra- 
tion of mm-sized chondrules in turbulent eddies followed by direct 
gravitational collapse ( Cuzzi et al.|2008 K the streaming instability 



i Youdin & Goodman 2005 i, and particle trapping in zonal flows 
Johansen et al. p009t 'r As demonstrated by [Johansen et al.|j2007| 
processes of this kind may ultimately lead to the rapid for- 
mation of bodies larger than Ceres via gravitational collapse of lo- 
cal concentrations of metre-sized objects, where the planetesimal 
masses depend on the disc mass. 

Owing to their size, the emergence of planetesimals marks the 
transition into the gravitationally-dominated stage of planet forma- 
tion. Runaway growth quickly leads to the formation of oligarchs 
( [Wetherill & Stewar t 1993; Ke nyon & Bromley|2009| , which sub- 
sequently grow into planetary embryos and cores via mutual colli- 
sions and accretion of smaller bodies jida & Makino|1993{|Kokubol 
|& Ida|199 8'). For this runaway growth to proceed efficiently, the ve- 
locity dispersion within the planetesimal swarm must remain sig- 
nificantly smaller than the surface-escape velocity of the accret- 
ing cores. Accounting only for self-stirring of the population, this 
is certainly expected in a fully laminar protostellar nebula. How- 
ever, the explanation of the typical accretion rates of T Tauri stars 
requires an anomalous source of viscous redistribution of angular 
momentum, most likely explained by disc turbulence. This poses 
the question whether or not runaway growth remains a realistic 
proposition in a turbulent protoplanetary disc. 

Planetesimals and protoplanetary embryos embedded in tur- 
bulent PPDs are subject to stochastic gravitational forces caused 
by density waves driven by the turbulent Maxwell stresses act- 
ing within the disc ^Nelson & Papa loizou"'20 04). Th is issue was 
addressed in our previous work ((Nelson & Gressel||2010| here- 
after "Paper I"), where we compared global and local magneto- 
hydrodynamic (MHD) simulations of protoplanetary discs involv- 
ing turbulence driven by the magneto-rotational instability (MRI, 
|Balbus & Hawley|1991^ . Following the evolution of swarms of em- 
bedded test particles, we studied the effects of magnetic turbulence 
on the dispersion of planetesimal semimajor axes and the growth of 
their internal velocity dispersion. In agreement with an earlier study 
( |Nelson|2005) , we found that fully-developed disc turbulence at a 
level consistent with observational constraints (see e.g. Armitage| 
|1998( and references therein) would cause large-scale diffusion of 
planetesimals over distances of several au (in contradiction of solar 
system constraints) and induce their collisional destruction. 

However, as originally contemplated by |Pneuman & Mitchell| 
( |1965^ , sustaining magnetic fields within the protosolar nebula re- 
quires the gas to be sufficiently ionised. This was later pointed out 
in the context of magnetised disc turbulence by |Safronov| ( |1972| l. 
In the wake of its discovery in the context of accretion discs, the 
linear development of the MRI was studied in detail for the case of 
ion-neutral drift (Blaes & Balbus 1994), and Ohmic diffusion (Jin 
|1996[[Sano & Miyama,1999) . The requirements for non-linear tur- 
bulence to be sustained have been studied by means of direct simu- 
lations for the case of Ohmic resistivity ( Fleming et al. 2000), and 
ambipolar diffusion ( [Hawley & Stone 1998, Bai & Stone 2011) . 
Because of the stabilising effect of epicyclic oscillations, a purely 
hydrodynamic flow is expected to be stable via the Rayleigh cri- 
terion ( |Rayleigh|I917| l. Lacking a plausible source for a sufficient 
number of free electrons, the ffow is likely to remain in a laminar 
state. 



Protoplanetary discs in the T Tauri stage are cold and dense, 
leading to the above mentioned low levels of ionisation ( |Ume-| 
ibayashr||l983^ . Based on the fact that the dominant ionisation 
sources are likely external to the disc itself, [Gammie| ( [T996) pro- 
posed a layered model for partially ionised accretion discs: While 
the disc surface layers are well ionised by stellar X-rays, and sus- 
tain a turbulent accretion flow caused by the MRI, the shielded disc 
interior is expected to harbour a 'dead zone' where the flow re- 
mains laminar implying negligible levels of accretion. Such lay- 
ered discs were first investigated numerically by [Fleming & Stone] 
(2003), who included a magnetic diffusivity varying with height. 
Their box simulations demonstrated that the flow within the dead 
zone-region maintains a modest Reynolds stress due to waves being 
excited by the active layers. 

Owing to the high gas density in the bulk of the disc and the 
embedded sub-micron dust grains, a number of physical effects 
may significantly affect the overall ionisation balance. This holds 
true for both PPDs in general, and the protosolar nebula in partic- 
ular ( |Hayashi|1981 1. The effect of dust on charge carriers was first 
investigated by Sano et al.|(2000). In a series of papers, [Tlgner &[ 
[Nelson[ studied the effect of small dust grains and different chemi- 
cal reaction net works jllgner & Nelson|2006a^ , and the role of tur- 
bulent mixing jllgner & Nelson[[2006b^ . An additional complica- 
tion arises when one considers the Hall term in Ohm's law (: Wardle[ 
[1999| >, introducing a dichotomy with respect to the mutual orienta- 
tion of any vertical magnetic ffux and the rotation axis ( jWardle &[ 
[Salmeron|2011| l. A study of the eflFect of the Hall term on the MRI 
turbulence saturation amplitude by [Sano & Stone[(2002 ^ found it to 
be little changed for the range of Hall parameters they investigated. 

Studies of the effects of ohmic dissipation determined by time- 
dependent reaction networks coupled to the MHD evolution have 
been presented by [Tumer et al.[ ( [2007) and ,Ilgner & Nelson,(2008) . 
These works have shown that enlivening the dead zone via turbulent 
mixing of charge carriers is only efficient in the absence of small 
dust grains. However, if a significant population of micron sized 
grains is present, the adsorption time scale of free electrons is short 
enough to maintain the dead zone ( [Turner & Sano[2008[[Bai[201 l| l. 
[Inutsuka & Sano[ ( [2005| l have suggested that fast moving electrons 
generated in magnetic reconnection events may provide a source of 
ionisation that can remove or modify the dead zone, but detailed 
calculations demonstrating the feasibility of this idea have yet to be 
done. 

Building on this line of work, we have extended our study 
from Paper I to the case of stratified PPDs. This includes a fiducial 
model being fully MRI active, and two models including a magnet- 
ically inactive midplane region ( [Gressel, Nelson & Tumer][2011| 
hereafter "Paper 11"). The key results of this study were that the 
stirring of the planetesimals' velocity dispersion was much reduced 
in the dead zone, as was their radial diffusion ~ satisfying the ob- 
servational constraint that radial mixing of asteroids has probably 
not occurred over distances ^ 0.5 au jGradie & Tedesco|1982[ l. 

This possibly allows for the continued growth of planetesi- 
mals rather than destruction or erosion during collisions, and it was 
suggested that dead zones may provide safe havens for planetesi- 
mals. One perceived limitation of Paper II was that we restricted 
ourselves to the case of a minimum mass protosolar nebula. While 
one would expect the strength of the stirring to scale approximately 
linearly with the disc mass (all other things being equal), the rela- 
tive level of density fluctuations might be affected by the extended 
width of the dead zone. It is the purpose of this paper, to eluci- 
date on the net effect of these two counteracting dependencies, and 



© 2002 RAS, MNRAS 000,[TpO| 



Disc mass and field dependence ofplanetesimal stirring 3 



Table 1. Overview of simulation parameters, including model Dl from Pa- 
per II. The fiducial model D 1 . 1 is identical to D 1 and serves as a control for 
the changes described in Sect.|2.3|below. 



10" 





Pmid 


cm- 


domain [//] 


resolution 


Dl/Dl.l 


1 


134.6 


3x12, ±5.500 


72x144x264 


D1.2 


2 


269.2 


3x12, ±5.667 


72x144x272 


DI.4 


4 


538.4 


3x12, ±5.833 


72x144x280 



moreover to examine the influence of varying the imposed mag- 
netic field strength. 

This paper is organised as follows: the physical model and nu- 
merical methods are briefly recapitulated in Section [2j for a more 
thorough description, we refer the reader to Paper II. Simulation re- 
sults are presented in three parts: in Section[3]we study the depen- 
dence on the disc mass/column density, and in Section|4] we show 
results on the variation with the external magnetic field permeating 
the disc vertically. A modified scaling relation to encompass all the 
obtained results is then developed in Section|5] Implications for the 
scope and applicability of our results for planet formation theory 
are discussed in Section|6] and conclusions are drawn in Section|7] 



2 METHODS 

The simulations presented in this paper are based on our fiducial 
dead zone model Dl discussed in Paper II. We again solve the 
standard equations of resistive MHD in a locally corotating, Carte- 
sian coordinate frame (x, y, z), adopting the shearing box formal- 
ism (see |Gressel & Ziegler|2007l and eqns. in sect. 2 of Paper II). 
Embedded in the dynamically evolving plasma are swarms of 25 
non-interacting massive test particles, which are affected by disc 
gravity and the inertial forces in the local frame of reference. All 
simulations utilise the mrvana-iii code ( |Ziegler|2004^ ; further mod- 
ifications to the solver are documented in sect. 2.4 of Paper II. 

The base state of our model is given by isothermal stratifica- 
tion in a fixed gravitational potential ^{z). Unlike in earlier studies 
( [Fleming & Stone||2003| |Oishi et al.][2557) , and to allow for suf- 
ficiently wide MRI-active regions, model Dl used a box size of 
5.5 scale heights, H, on each side of the disc midplaneQSased on 
our previous analysis of the influence of box size on gravitational 
stirring (Paper I), we furthermore maintain a horizontal box size 
of 3 X 12 H, permitting the excitation of spiral density (SD) waves 
( |Heinemann & Papaloizou|2009a|b||20T7| . 

Because the resistivity r] depends on a physically motivated 
ionisation model involving chemical rate equations, we have to in- 
troduce a unit system to convert code units into standard units. To 
maintain continuity with Paper II, we adopt the same basic disc 
model, which is similar to the widely used "minimum-mass" pro- 
tosolar nebula ( |Hayashi|1981^ . At 5 au the Hayashi model has a col- 
umn density Z = ISOgcm""^ and a disc opening angle Hjr = 0.047. 
We adopt similar values E = 135 gcm"^ and H/r = 0.05, giv- 
ing a temperature T = 108 K, and an isothermal sound speed of 
c's = 667 ms"' (also cf. Tab. [5}. Because the minimum-mass as- 
sumption does in fact only specify a lower limit to the expected 



With our definition of H, the equilibrium den sity profile is proportional to 
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Figure 1. Vertical profiles for the runs presented in this paper. We plot 
the magnetic Reynolds number, Rm = C1H~ 7j(zT^ determined from our 
ionisation model. The three line represent the initial models for runs Dl, 
D1.2, and D1.4, respectively (cf. Tab.[lJ; short vertical bars illustrate the 
varying domain size L- . 



Table 2. List of model parameters. Note also that magnetic field values for 
each run are listed in the third column of table|3] 



parameter 


symbol 


value 


unit 


orbital location 


a 


5.0 


au 


disc aspect ratio 


h 


0.05 




local temperature 


T 


108 


K 


mean molecular weight 




2.0 


amu 


gas/dust ratio 




0.001 




dust particle size 




0.1 


lira 


dust grain density 


Pi 


3.0 


gcm"^ 


disc column density 


I 


135 -538 


gcm"2 


magnetic field strength 


6; 


2.68 - 43 


mG 



exp(-z /(2//-)), i.e. H is identical to the 'h' in Okuzumi & Hirose|^201 1^ 



mass of the protosolar nebula at the age when planetesimal forma- 
tion likely occurred, we here study further models with double and 
four- times the midplane density pmid (see Tab. [TJ. The ionisation 
fluxes for X-rays (XR), short-lived radionuclides (SR), and cosmic 
rays (CR) are unchanged from model Dl, i.e., XR and CR fluxes 
are the nominal values, and SRs are enhanced by a factor of ten; see 
Sect. |2.2.T| for further details. For studying the influence of the net 
vertical field (NVF), we performed additional runs Dl-WF (with a 
weak field), Dl-NVF a/b (with varying field), and two extra models 
based on D 1.4 with weaker fields. We adopt a standard resolution 
of 24 grid cells per H in the radial, and vertical directions; for the 
azimuthal direction we use a reduced grid spacing of 12/ H. A basic 
convergence check is provided in Appendix [A] 

2.1 Initial and boundary conditions 

A primary aim of this work is to construct discs of varying mass 
in which the active zone remains essentially unchanged in its mass 
and level of turbulence in each model, with the only change being 
in the mass and vertical thickness of the dead zone. This component 
of our study then becomes one in which we examine the effect of 
varying the mass of the dead zone keeping all other parameters 
fixed (see Tab.|2]for reference). The numerical set-up we adopt for 
the models of varying disc mass is designed to achieve this aim. 

The initial magnetic field is a superposition of two com- 
ponents, a standard zero-net-flux (ZNF) configuration S-(x) ~ 
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sin(2;?r x/L^), and a weak additional vertical net-flux (NVF) compo- 
nent representing an external or dynamo-generated field. We vary 
the strength of the initial ZNF flux to produce a plasma parameter 
(i.e. the ratio of thermal to magnetic pressure) /3j, - 50 independent 
of the midplane density (and hence pressure) of the three studied 
models. In contrast, we keep the absolute strength of the external 
net vertical flux (which is conserved during the simulations) fixed 
at a value corresponding to 60 = 10.73 mG in physical units, result- 
ing in a weaker relative field in the midplane of models with higher 
mass loading. Note that the NVF component leads to a plasma pa- 
rameter varying with height. As a consequence, the region of inter- 
mediate field strength, i.e. the vertical range, where the most unsta- 
ble nonaxisymmetric MRI modes fit into the box, moves up in the 
disc. We expect the resulting active regions to be of similar vertical 
extent irrespective of the mass loading, but separated by a wider 
dead zone in the case of heavier discs as is illustrated by shaded 
areas in Fig. [T| To maintain near-identical active zones across all 
runs, we adjust the vertical box size accordingly (cf. Tab.[TJ. 

The simulations that examine the influence of varying the 
strength of the net vertical magnetic field strength (models Dl- 
NVFa and Dl-NVFl)) are initiated with specific values of the mag- 
netic field (Bq = 2.68 mG and Sq = 10.73 mG, respectively), and 
after every 40 orbits the field strength is incremented by factors of 
two up to some maximum value. Although we only adopt two run 
labels for these simulations, these simulations in fact explore the 
influence of six different external magnetic field strengths ranging 
between 2.68 mG and 86 mG. 

As discussed in more detail in Paper II, for the fluid vari- 
ables, we use vertical boundary conditions allowing material to 
leave the box but prevent inflow. Mass loss associated with the out- 
flow boundaries is compensated by continually re-instating the ini- 
tial density profile in each grid cell by means of an artificial mass 
source term (cf. [Hanasz et al.||2009) . Because of the low density 
in the halo, the relative change in the disc mass due to outflow is 
very low (on the order of 10"* per orbit). We have verified that the 
additional term does not affect the obtained stirring amplitudes. As 
in our previous work, we compute a locally and temporally vary- 
ing magnetic diffusivity t](x, t). This adds to the level of realism 
achieved in earlier simulations of dead zones by [Fleming & Stone] 
(j2003j and |Oishi et al.| ( [2007^ , who used static diffusivity profiles. 
Here we briefly re-capitulate the physical ionisation model we im- 
plement and refer the reader to section 2.3 of Paper II for a more 
detailed discussion. 



2.2 The difTusivity model 



Because of the expected dominance of small dust grains jSano et aT] 
[2000 ; Ilgner & Nelson 2006a), we decide to avoid following the de- 
tailed non-equilibrium chemistry and adopt a simplified approach 
for the gas-phase reactions. Assuming that recombination happens 
much faster than any dynamical mixing timescale in our system, 
we update 77 according to a precomputed table derived from the re- 
action network in inodel4 of |Ilgner & Nel son (2006a). When com- 
puting the resistivity, we include the contributions of all the charged 
species following |Wardle| ( |2007^ , eqns. (21)-(31). For the free pa- 
rameters, we assume the same values described in Paper II, i.e. dust 
grains of size 0.1/jm and with density 3gcm"^ and a dust-to-gas 
mass ratio of 10"^ (such that we tacitly assume that 90 percent 
of the solid grains have already grown to become larger bodies). 
The gas-phase abundance of Magnesium is taken to be depleted 
by a factor lO"** compared to its solar value (with the remainder 
assumed to be bound-up in grains). The key factor governing the 



diffusivity 77, is the local ionisation rate f (x, t), which is computed 
from the external irradiation by evaluating column densities to both 
the upper and lower disc surfaces. Because of the expected radial 
density features, related to the excited spiral waves, this is done on 
a per-grid-cell basis. 

2.2.1 Ionisation sources 



The ionisation model is founded on the work of [Turner & Drake! 
l |2009j l, who studied how stellar X-rays, radionuclides, and ener- 
getic protons can influence the shape and extent of a possible dead 
zone via their effect on Ohmic diffusion. Implementing the model 
of |Turner & Drake| ( (2009l l, we focus on stellar X-rays (XR), and in- 
terstellar cosmic rays (CRs) as the prime sources of ionisation. As 
a representative value, we chose Lxr -2x10'" erg s"' as suggested 
by Garmire et aL]|2000| ). Applying a simple fit to the Monte-Carlo 
radiative transfer calculations of |Igea & Glassgoid|jl999f , we ap- 
proximate the ionisation rate due to X-rays by 



(TxR = 2.6xl0"'^s-' e"^"'^™ ■ 



(1) 



with Ea and the gas column densities above and below a given 
point, and Sxr = 8.0gcm~^ the assumed X-ray absorption depth. 
For the vertical attenuation of interstellar CRs illuminating the disc 
surfaces, we adopt the formula given in |Umebayashi & Nakano| 
l |2009) : 



1 + 



(2) 



where Scr = 96 g cm is the cosmic ray attenuation depth 1 Ume 



[bayashi & Nakano|I981[ l, and dots indicate the corresponding con- 
tribution from the second column density S;,. We further include 
an ambient ionisation due to the decay of short-lived radionuclides 
(SR). As already done for model Dl, and to somewhat limit the dy- 
namic range in 77, the related ionisation is chosen lOx the nominal 



value of fsR = 3.7x10"'^ s"' quoted in Turner & Drake 1 2009 1 



2.3 The induction equation 

When simulating protoplanetary discs with dead zones using the 
shearing box approximation, there are two issues that arise which 
must be addressed to increase the feasibility and realism of the 
models. The first occurs because we are solving an advection- 
diffusion equation in time-explicit fashion. The time step size is re- 
stricted by the time for diffusion across one grid cell, and because 
this quantity scales with the square of the grid spacing, the pres- 
ence of large diffusion coefficients can render the computational 
effort prohibitive in highly resolved runs. This was addressed in our 
previous paper by limiting the dynamic range in r] and applying the 
technique of temporal sub-cycling to the diffusive part of the induc- 
tion equation (cf. sect. 2.3.2 in Paper II). With the more extended 
dead zones of models D1.2 and D1.4, the time scales become even 
more disparate, and we have developed a modified scheme to deal 
with this, described in detail in Appendix [B] 

The second issue relates to the generation of strong azimuthal 
magnetic fields through the winding up of net radial fields in the 
box. Although our simulations begin with zero net radial magnetic 
field in the computational domain, advection of horizontal field 
components through the open vertical boundary leads to the gener- 
ation of net radial fields that may diffuse into the dead zone. Once 
there, Keplerian shear may wind them up to generate strong az- 
imuthal fields that subsequently leak back into the active regions. 
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modifying the turbulence there. This issue has been discussed by 
[Turner & S ano ( 2008 ), where it was noted that the winding of fields 
generates a strong radial gradient in the azimuthal field strength due 
to the Keplerian rotation profile in global discs, leading to radial 
diffusion of the strong fields that limits their growth. This radial 
diffusion cannot occur in a shearing box due to the uniform shear 
and periodic boundary conditions (but has been observed in the 
global simulations presented by |Dzyurkevich et al.| l j2010 )), so we 
have developed an approximate scheme designed to crudely mimic 
the radial diffusion of fields in global discs. We refer to this scheme 
as super-box scale diffusion, and we describe its implementation in 
Appendix [B| 



3 DEPENDENCE ON DISC MASS 

In Paper II, we restricted our analysis to the fiducial case of a min- 
imum mass protosolar nebula, and this naturally begs the question 
how the results will be affected if one looks at more massive discs. 
A major aim of this study is to compute disc models with differ- 
ent masses/surface densities that have very similar active layers in 
terms of mass and turbulent activity. The main diff'erence between 
the models will then be the physical depth and mass of the dead 
zone. Density waves that are excited in the active layer will propa- 
gate into the dead zone and induce gravitational stirring of the plan- 
etesimals there through the density fluctuations that they generate. 
A key question is how this stirring depends on disc mass. 



3.1 Scaling relations for linear density waves 

Based on the results of simulations presented in paper II, we expect 
the disc models described below to have the following midplane- 
symmetric three-layer structure: a non-turbulent dead zone located 
at and above the midplane in which the MRI is quenched by Ohmic 
resistivity, on top of which lies a turbulent region in which the MRI 
operates unaffected by resistivity. At disc altitudes above the MRI- 
active region lies a magnetically dominated halo in which the field 
strength is too large for the MRI to operate because unstable modes 
do not fit within the available vertical extent. Density waves excited 
at the interface between the dead zone and MRI-active layer will 
propagate into the dead zone, creating density fluctuations at the 
midplane. 

We denote the density at the active/dead zone interface as 
Pint- The disc models of different mass we compute are designed 
to generate discs in which very similar active layers lie above 
and below dead zones of different mass, and for which the mid- 
plane density, p,„id, scales with the disc mass. We expect each 
of our models therefore to have very similar values of the den- 
sity. Pint, and perturbed velocity, (5vi„,, at the dead/active zone in- 
terface. Linear density waves that travel sonically into the dead 
zone should conserve energy as they propagate, such that we ex- 
pect Pn,id(^^'mid)^ - Pin((<5vint)^> where (5v„iid is the perturbed velocity 
at the midplane. Relative density fluctuations for isothermal sound 
waves satisfy the relation 



Sp 
P 



6v 



(3) 



such that relative density fluctuations at the disc midplane should 
obey the scaling 



Translated into absolute density fluctuations, this implies a scaling 
^Pmid ^ yjPmii- H is reasonable to expect the fluctuating torque, 
Fy, induced by density fluctuations near the midplane will scale 
linearly with 6p,^\i, leading to the naive expectation that stochastic 
torques will scale according to F,, oc y[p^. 

At the time of writing an extensive study of MRI turbulence 
in layered accretion discs appeared in print ( [Okuzumi & Hiro^ 
|201 \) . For a broader discussion about scaling relations in disc mod- 
els with dead zones we refer the reader to this recent publication. 



3.2 Disc evolution 

As in sect. 3.2 of Paper II, we will start our discussion by looking 
at the evolution of the gas disc as shown in Figure [2] Owing to 
the additional magnetic diffusion term (cf. Sect. |B2| irthe overall 
evolution of model Dl.l (left hand panels) appears somewhat less 
intermittent than the identical setup from model Dl in Paper II. 
This is also true for the heavier models D1.2 (centre) and D1.4 
(right hand panels). 

The different models in fact appear very similar, e.g., in terms 
of the cycle frequency of the field reversals and the overall field am- 
plitude. Also the stresses and relative density fluctuations within the 
MRI-active regions are nearly identical. The only apparent trend 
is seen in the relative density fluctuations in the midplane region, 
which fall off with mass loading. The relative amplitudes, as listed 
in column 7 of Tab. |3] are 0.091, 0.058, and 0.038, respectively, 
which implies a scaling with the midplane density to the power 
of ~ -0.63 ± 0.01; this is in decent agreement with the expected 
scaling ((5p„,id/Praid) °^ Pm!d^ discussed above in Sect. 3.1 The dis- 
crepancy can partly be ascribed to the varying 6p |a=i at the layer 
interface (cf. Sect. |3.3| below), which shows a slight trend ~ with 
a power ~ -0.05 ± 0.01 - towards weaker absolute fluctuations 



I P /mid 
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in heavier discs. Overall, the results indicate that the propagation 
of density waves from the transition region between the active and 
dead zones occurs without significant dissipation of energy. 

Because our model assumes that the main source of ionisation 
is exterior to the disc surface, the width of the dead zone depends 
on the column density of the shielding material. The extent of the 
dead zone as a function of various ionisation sources was studied 
in detail by [Turner & DrakeH2009j . In our simulations, we deter- 
mine the transition region via the criterion that the Elsasseij^Jnum- 
ber A = v^KClr]) equals unity (indicated by a solid black line in 
Fig. |2j. This yields a time-averaged half thickness of 1.66±0.04, 
2.02±0.03, and 2.33±0.03// for models Dl.1-4, respectively. In 
other words, doubling the disc mass increases the extent of the dead 
zone by a third of a pressure scale height H. As can be seen from 
Fig. |2] the extent of the different layers is remarkably stationary. 
As a result, the dead zone thickness only fluctuates on the per-cent 
level. This has, of course, to be taken with a grain of salt. Exter- 
nal conditions (e.g., gas inflow from the parent molecular cloud, 
stellar X-ray luminosity, intensity of arriving CRs, . . .) in a realis- 
tic protoplanetary disc are likely to vary strongly both in time and 
space. However, the relative quiescence of our experimental setup 
will allow us to quantify the secular evolution of the embedded 
planetesimals more accurately. 



^ We note, in passing, that this definition is formally identical to what has 
been called a Lundquist number (or even magnetic Reynolds number) in 
related publications. Irrespective of its name, A < 1 serves as a robust indi- 
cator for suppression of the fastest-growing MRI mode (see |Pessah|2"0l0) . 
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Figure 2. Space-time evolution of the mean toroidal field (upper row), total Maxwell stress M,^ = - (Br B^) (middle row), and relative density fluctuation 
dp/p, (bottom row). The different panels show runs Dl.l (left), D1.2 (centre), and D1.4 (right). The interface between the live and dead zones is indicated by 
a solid line representing the criterion A = v^j^/(Slri) = 1 . Note the trends of a wider dead zone, and weaker relative density perturbations when going to higher 
densities, while most disc properties remain virtually unchanged. 



3.3 Disc vertical structure 

Having assured ourselves of the stationary nature of the flow, we 
can now proceed and look at the characteristic structure of the lay- 
ered disc. In Figure [3] we plot the time-averaged (20-220 orbits) 
Maxwell stress M^^ = -{Bi-B^) for the three models D 1.1-4, of 
varying column density. As discussed in the last section, owing to 
the random diffusion of flux into the midplane region, there is no 
systematic to the residual level within the dead zone. As noted ear- 
lier, we expect the MRI-active zone to move away from the mid- 
plane for the higher density models, and this trend is clearly seen. 
As a result, we infer very similar vertically-integrated mass accre- 
tion rates of 7.90 x 10"**, 7.97 x 10-^ and 7.86 x lO"** yi-' for 
runs Dl.l -4, respectively. This is in-line with the volume-averaged 
values stated in column 4 of Tab. [5] and demonstrates that higher 
mass loading does not affect the level of activity or mass in the ac- 
tive zone, but merely the vertical extent and mass of the dead zone, 
such that the nonlinear evolution of the disc models arising from 
the specified initial conditions is as intended. 

Figure |4] shows the same shift away from the midplane for 
higher disc mass in the hydrodynamic stresses (pVrSv^). The as- 
sociated, vertically integrated mass accretion rates are 1.28 x 10"*, 
1.33 X 10"^ and 1.41 x lO^^Moyr"' for models Dl.1-4, respec- 
tively, i.e. the bulk of the transport is due to magnetic stresses. 
Within the dead zone, the Reynolds stresses scale with the gas den- 
sity but this trend is very weak at their respective peak position, 
which lies just outside the layer interface. This region is of partic- 
ular interest as it marks the position in the disc where the turbulent 
fluctuations created by the MRI-active layers are most energetic in 
terms of their mass loading and hence their associated momenta. 
Tracing the interface via the A = I criterion, we obtain very simi- 
lar absolute density fluctuations at this line. The respective values 
are 6p |a=i = 0.045, 0.043, and 0.042 for models Dl.I-4 

We conclude that, in all three models, the dead zone sees a 
very similar boundary separating it from the active zone, with the 
only difference being in its physical separation from the particle 
swarm located near z = 0. In the following section, we will deter- 
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Figure 3. Time-integrated Maxwell stresses for models Dl.1-4 as labelled. 
The vertical displacement of the curves demonstrates nicely why all models 
produce the same mass accretion rate. Shaded areas indicate the extent of 
the dead zone in each model. 




c 10 



, y 

, , I , , , i1 o 

-2 2 
height z [H] 

Figure 4. Time-integrated Reynolds stresses for models Dl.1-4. The same 
vertical displacement as in Fig. [3] is seen. The peaks in the stress sits just 
outside the transition region between the active and inactive layers (shaded). 
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Figure 5. Torque distributions F,. for tlie models Dl.1-4 witli varying disc 
mass. Residuals indicate tlie deviation from a normal distribution. 

mine the consequences this has on the gravitational forcing exerted 
onto the planetesimals, and examine if the stochastic torques follow 
the scaling predicted in Sect. |3.1[ 

3.4 Gravitational torques 

It was established in the global simulations of |Nelson & Papaloizou| 
l |2004^ and |Ne lson ( 2005 1 that density fluctuations from developed 
turbulence lead to stochastic gravitational forces which can have a 
significant impact on the dynamical evolution of embedded plan- 
etesimals and protoplanets. We remark that earlier shearing box 
simulations of MHD turbulence have shown the development of ax- 
isymmetric pressure bumps or 'zonal flows' ( [Johansen et al.|2009| l, 
which may influence the stirring of planetesimals. Analysis of our 
simulation D 1.1 suggests that the midplane density supports a long- 
lived (20-30 orbits), axisymmetric density asymmetry across the 
radial width of the box with amplitude ~ 1%. This may be an in- 
dication that weak zonal flows arise even in the presence of a dead 
zone. The role of this in stirring planetesimals is unclear at present 
and will the subject of a future study. 

In our previous papers, we have quantified the level of stirring 
for a range of disc models. We here perform a similar analysis and 
find broad agreement with previous runs. In particular, model Dl.l 
agrees well with the earlier run Dl from Paper II - we attribute the 
discrepancy of about ten per-cent in the torques between Dl (0.06) 
and DI.I (0.054) to the lower level of sampling noise achieved in 
the new set of simulations, where a more conservative choice for 
the gravitational softening parameter was taken. We note that when 
statisticalljl^correcting the original results from Dl for a weak sys- 
tematic dependence on the particle position on the grid, excellent 
agreement can be obtained. 

3. 4. 1 Torque distribution functions 

Torque time series, Ty{t), are derived from the gravitational forces 
acting on the particles along the y direction within our simulation 
box, and are recorded along the evolving trajectory of the parti- 
cles. Torque amplitudes, cr(ry), for the different runs are obtained 
from a Gaussian fit to the histograms as depicted in Figure [5] Av- 
eraging over the swarm of 25 particles and sub-intervals in time 

^ As different locations ought to be indiscriminate in a stochastic sense, this 
is done by binning the full data with respect to particle position (at sub-grid 
resolution) and then removing the median torque for each bin. 
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Figure 6. Autocorrelation functions (ACFs) of the torque time series, aver- 
aged over 25 particles, and 4 sub-intervals in time (shaded areas). The first 
and second zero crossing are indicated by a triangle. Labels give the fitted 
coherence time along with the relative amplitude of the wavelike feature. 

allows us to estimate error bars, and we infer values of aiT^) = 
0.054±0.004 for model Dl.l, and 0.052±0.003 for both models 
D1.2 and D1.4 (these values are also listed for reference in column 
eight of Table |3j. This demonstrates that the width of the torque 
distribution is remarkably insensitive to the column density, and 
models Dl.l through D1.4 agree to within five per-cent, which is 
somewhat smaller than the error bars. This result is in clear conflict 
with the simple scaling law based on conservation of wave energy 
for isothermal sounds waves (see Sect. |3.1^ . 

Naive expectations suggest that fluctuating torques experi- 
enced by solid bodies located near the disc midplane will scale 
linearly with absolute density fluctuations there, where we both 
predict and observe the scaling Spmii °^ s/Pmid- Thus, we expect 
the fluctuating torque amplitude to scale with the square-root of 
the disc mass in these simulations. As discussed in Sect. |3.3| we 
find that all models agree on the relative density fluctuations at the 
transition between the dead and active zone, and the absolute and 
relative density fluctuations at the midplane are in good agreement 
with the scaling arguments (also cf. Fig. [T3j. This suggests that a 
more subtle effect is responsible for modifying the expected linear 
scaling between the midplane density fluctuations, Spn^id, and the 
rms torque fluctuations, (T(r,,), shown in Figure [5] A factor of four 
change in disc mass between models Dl.l and D1.4 should lead 
to a factor of two change in the torque amplitude, and this would 
clearly be detectable if present. We defer a more detailed discussion 
of this issue until Sect.|5] where we provide evidence to resolve the 
discrepancy between the simulation results and linear scaling. 

3.4.2 Torque autocorrelation 

How does the modified amplitude of the density waves influence 
the temporal coherence of the stochastic torques? We recall from 
our previous work that the torque autocorrelation function (ACF) 
appeared as a superposition of a truly stochastic part and a mod- 
ulated "wavelike" feature (due to the quasi-periodic passage of 
density waves through the particle's location). This dual character, 
which is also seen in the corresponding figs. 9 and 15 of |Oishi et al.| 
(j2007j and |Yang et al.|(2009j ), was described in terms of a formula 
which we introduced in section 3.3 of Paper I: 

5r(r) = [(l-a) + a cos(2n cot)] a''^' . (5) 

Fitting parameters are the relative amplitude a, the frequency of the 
wavelike modulation u>, and the coherence time of the stochas- 
tic envelope. The resulting fits for the models Dl.1-4 are plotted in 
Fig.|6] where we also state the obtained values for the relevant fit- 
ting parameters. In terms of the parameter of primary interest, the 
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Figure 7. Random-walk behaviour for the rms radial displacement <t(/S.x). 
The right hand axis shows the absolute dispersion o"(Aa) at ao = 5 au. Note 
that model D1.4 appears to show a slightly reduced tendency for particle 
dispersion. 



coherence time t^, models Dl.l and D1.2 show excellent agree- 
ment with model Dl from Paper II. Note that while model Dl 
showed a ~ 55% wavelike modulation, this is now increased to 
~ 75%. We believe that this relative increase in the wave feature 
is, in fact, a reduction of the stochastic component caused by the 
lower level of random sampling noise in the new set of simulations. 
When correcting for this in model Dl, we also obtain a wavelike 
amplitude of ~ 75%. The modulation is even more pronounced in 
model D1.4, which moreover shows a significantly reduced coher- 
ence time of Tc = 0.22. One can see by inspection of the over- 
plotted model curves that there exists some tension between our 
simplistic approach and the actual data. We conclude that a more 
accurate determination will require a more thorough understanding 
of the ACF shape. Yet, ultimately the derived stochastic description 
via criXy) and can be verified by comparison with the amount of 
particle diffusion present in the simulations, as we shall demon- 
strate in the next section. 

3.5 Radial diffusion of planetesimals 

Having demonstrated that the gravitational torques acting on em- 
bedded planetesimals can be approximated by a normal distribu- 
tion, and their temporal correlations possess a finite coherence time, 
we can now exploit these two properties to predict the amount of 
particle diffusion based on a Fokker-Planck-type model. This ap- 
proach was described in detail in section 6 of Paper II, and we 
here briefly recapitulate the essential ingredients. Particle spread- 
ing occurs because stochastic torques cause diffusion of planetes- 
imal angular momenta at a rate given by the diffusion coefficient 
£)p = criTyf- T^. This leads to an estimated spread in semi-major 
axis a after an evolution time, t, of 



Afl ■■ 



(6) 



as was derived in eqns. (16)-(19) in Paper II. This estimate for Aa 
can then be confronted with a random walk-type fit to the rms par- 
ticle displacement in semimajor axis due to diffusion, a(t^x), ob- 
tained in the simulations 



cr(Ax) 
H 



■ C^(Ax)Vf, 



(7) 



corresponding to eqn. (20) in Paper II. The rms spread in the local 
displacement = A(a - ag) is plotted in Fig. |7] for the differ- 
ent models. Owing to the relatively small number of particles, the 



curves display substantial random fluctuations, yet the general be- 
haviour of a random walk-like growth emerges, and it can further 
be seen that the rate of planetesimal spreading in the runs is sim- 
ilar, as expected from the earlier discussion concerning the simi- 
larity of the stochastic torques in the three runs. The fitted values 
obtained for C^(Ax) are (4.56 ± 0.54) x lO"'' for model Dl.l, and 
(4.36 ± 0.32) X 10"* for model D1.2, in good agreement with the 
value for run Dl from Paper II, which is listed in Tab. [3] As can 
also be seen from Fig. [t] the value of (2.98 ± 0.77) x 10"" for model 
D1.4 is somewhat reduced. This is consistent, however, with the 
lower inferred coherence time in this case. 

To be more specific, let us compare the prediction based on |6| 
with the measurement of Co-(Ax), which can be related via |7j. For 
model Dl.l, we infer a spread of Aa = 0.00161 au at a distance of 
ciQ = 5 au, and after a run-time of t - to =200 orbits. This value 
may readily be checked by reading it off the right hand axis of Fig- 
ure |7] By means of specifying t^. and cr(rv), the diffusion equation 
predicts Aa = 0.00148 au, which is only about ten percent smaller 
than the actual value. A similar comparison for model D1.2 yields 
Aa = 0.00154 au from Fig. |7] opposed to an estimated 0.00140 au 
from the diffusion approach, i.e., agreement again to within ten per- 
cent. Finally, for model D1.4 we predict Aa = 0.00122 au via (j6]l, 
and observe a value of 0.00105 au, which is about 15% smaller than 
the estimated value. Compared to the results from Paper II, this 
means a substantially improved accuracy of the predictions. We at- 
tribute this to the more stationary behaviour of the new simulations 
effected by the inclusion of super-box scale dissipative effects (cf. 
Sect. |B2[ ). We conclude that the simplified description of particle 
spreading in terms of a diffusion equation can be used as a reliable 
tool in predicting the level of dispersion in the immersed planetesi- 
mal population within a dead zone. Moreover, the good agreement 
between the estimated and measured values lends support to the 
accuracy of the coherence times fitted via (jsj. 

3.6 Eccentricity stirring 

To complete our discussion on how the embedded planetesimals are 
affected by the turbulence, we now briefly consider the excitation of 
their eccentricity. This has been done in some detail in section 5 of 
Paper II, including a discussion of the long-term evolution and re- 
lated saturation mechanisms. In this section, we consider the results 
of the simulations directly, and leave discussion about their impli- 
cations for planetesimal accretion and planet formation theory until 
later in Sect. [6] of this paper. 

Figure [sjshows that the rms eccentricity curves for the differ- 
ent models resemble a braided band. Owing to the modest num- 
ber of embedded particles, the chance fluctuations are quite pro- 
nounced again. We register an overall spread of the curves of about 
a factor of two, which explains the scatter in the fitted amplitudes, 
Co-(e), describing the random-walk curve 



(r(e) r 
--}-^=CAe)<t. 
H/r 



(8) 



The values are listed in the last column of Tab. [3] and as a whole 
agree markedly well with the result from the earlier run Dl from 
Paper II (see Table 3 in that paper). The striking agreement of the 
three very different simulations Dl.l -4 in terms of their effect on 
embedded particles again highlights the rather stark disagreement 
with a scaling proportional to the midplane density fluctuations 6p. 
This will be discussed in detail in Section [5] and the implications 
of the obtained stirring relations for planetesimal evolution will be 
discussed in Sect.|6] Having established the dependence on the disc 
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Figure 8. Comparison of the rms random-walk eccentricity cr(e), of .sets 
of 25 particles in runs Dl.1-4. The dashed line shows the expected 
behaviour. The second ordinate indicates the radial velocity dispersion 
amongst the particles in absolute units. The stirring amplitude is virtually 
identical in all three runs. 



mass, we now move on to the effect of the external net flux perme- 
ating our simulation box. 



4 DEPENDENCE ON EXTERNAL MAGNETIC FIELD 

In the previous discussion, as well as in Papers I and II, we have 
tacitly assumed a fiducial value of about lOmG for the net verti- 
cal magnetic field permeating the accretion disc (see third column 
in Tab.[3j. This particular value was chosen to produce a turbulent 
transport coefficient compatible with commonly accepted T Tauri 
accretion rates. In a global picture, such a field naturally arises from 
local fluctuations - i.e., it would be hard to imagine that all the 
various subsections of a disc manage to retain zero vertical flux 
at the same time. [Sorathia et ari ( |2010[ l have recently verified this 
notion by dissecting global zero-net-flux simulations into shearing- 
box size portions, and found that the turbulent stresses from the 
emerging dynamo field match the scaling derived from local box 
simulations with an equivalent imposed vertical net-flux. While 
their simulations did not include explicit dissipation terms, there is 
no obvious reason why this correspondence should be lost in global 
models including dead zones jDzyurkevich et al.|2010^. Note , how- 
ever, that from unstratified local simulations Fleming et al.|pOOO) 
found a critical Rmc - Iff* necessary to sustain MRI in a zero- 
net-flux configuration. This threshold was recently confirmed for 
stratified discs by |Oishi & Mac LowH201 ll . Looking at Fig.[T] we 
see that this criterion is only met above about 3.5 // in our models, 
leaving very little space for a self-sustained dynamo mechanism. 
Because convergence is much harder to obtain without a net flux 
(see e.g. [Davis et al.pOlO] and references therein), we consider 
only runs that include net flux vertical fields in this work. 

As an alternative to a dynamo-generated disc field, one might 
envisage an external stellar field permeating the disc. Given the 
highly dynamical star-disc interaction during the T Tauri phase, 
strong magnetic fields are to be expected, and observations find 
stellar surface fields of roughly 2 - 3kG jBouvier et al.|[2007[ l. 
Assuming a dipolar field geometry implies a scaling of S = 
B^. (R^,/ry, and taking B^, = 3kG, and R^, = 2i?e, we arrive at 
~ 2.4 mG, and ~ 20yuG at r = 1 au, and r = 5 au, respectively. 
Clearly, the contribution of the central star is very moderate at 5 au, 
but should be considered when looking at models further in. In the 
following, we will ignore the actual origin of the precise value of 



the imposed field, and simply study how a varying level of the field 
strength influences our results. 



4.1 Preliminary considerations 

On theoretical grounds there are two effects which determine the 
expected level of turbulence: (i) the unmistakable scaling of the 
turbulent stresses with the net flux, and (ii) the extent to which tur- 
bulent activity in the MRI-sustaining regions can penetrate into the 
poorly ionised layer near the midplane. As for the first dependence, 
it is now safely established that the saturated state of the MRI de- 
pends critically on the net field ( [Hawley et al.|1995{ [Pessah et al.| 
|2007| l, and as [Yang et al.|p009l l have shown from unstratified sim- 
ulations without dead zones, this should translate directly into the 
level of turbulent stirring (cf. their fig. 3). In the following, we will 
set out to explore how this scaling is affected by a dead zone and 
try to establish the limits of its applicability. 

In the case of a net vertical flux (NVF), the most powerful 
MRI modes are so-called channel flows, i.e. axisymmetric solutions 
of the perturbed magnetic and velocity fields, only depending on 
the vertical coordinate z- It is a remarkable property of these modes 
that they remain exact solutions of the underlying equations far into 
the non-linear regime; this property even persists in stratified discs 
( [Latter et al.|2010j . Given the large-scale nature of these dominant 
MRI modes, the question arises whether the flow pattern can effi- 
ciently be broken-down into unordered motion. It is commonly ac- 
cepted that this transition into turbulence happens via parasitic in- 
stabilities (Goodman & Xu 1994; Pessah & Goodman 2009; Latter] 
[et al.[2009) feeding on the unstable MRI modes. A detailed analysis 
of how these modes look in the case of a stratified disc can be found 
in sect. 4 of |Latter et al.[j20l"0) . The key findings of this analysis are 
that (i) the parasitic modes have wavelengths of about /iMRi/2, and 
(ii) their growth rates are considerably reduced compared to the un- 
stratified case. The latter finding leads to the conclusion that chan- 
nel flows are expected to be persistent, if not dominant, features in 
stratified discs. This notion is supported by the 'streaks' seen in the 
density fluctuations (lower-most panels in Figs.[10[and[ll|(, which 
are remarkably similar to the pinching effect in fig. 3 of Lat ter et al.[ 
( [2010^ . Further evidence for the strong presence of channel modes 
comes from the horizontally averaged radial and azimuthal veloci- 
ties ii,-, and (not shown), which exhibit a remarkable correlation 
with S, , and B^, and notably show the opposite parity compared to 
their magnetic counterparts. This symmetry is highly indicative of 
the solutions derived in [Latter et al.[ ( [2010^ , cf. their fig. 1, and a 
direct consequence of the underlying set of equations. To conclude 
this digression, we note that while the presence of MRI channels 
was already mentioned in Paper II, their appearance is even more 
pronounced when going to higher field strength. 



4.2 Turbulence amplitude as a function of net-flux 

In the left panel of Fig.|9] we plot the time- and volume-averaged 
transport coefficients obtained from runs Dl-NVFa/b, covering a 
range of B, = 2.68 - 43 mG (i.e. field strengths spanning a factor 
of sixteen). We want to point out that this interval roughly coincides 
with the practical limits set by the MRI wavelength. 



:l(z) : 




(9) 



of the fastest-growing linear mode. At the weak-field end, the wave- 
length is just long enough to be well-resolved, while at the strong- 
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Table 3. Overview of simulation results including runs from papers 1 & II. For reference, we list the net vertical field ( ) (col. 3) which serves as an input 
parameter. All turbulent stress parameters (cols. 4-6) are normalised by the midplane pressure of model Dl, i.e., po = 1- Relative fluctuations Spip are midplane 
values. The width tT{Yy) of the torque amplitude is obtained from Fig.|5] correlation times are from the ACFs in Fig.|6] and Ca-((^x) and Co-(e) are derived 
from Figs.|7]and[8] respectively. Note that runs NVFa/b apply dift'erent Tjiff (cf. Eqn. jB3j in Sect.|B2j, assuming = njr, and Injr, respectively. 
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Figure 9. Left: Scaling of the dimensionless transport coefficients T,-^ = 
(pv,-6v^ ), and M,-^ = (~ SB^SB^ ). The stresses roughly scale with B- to 
the power of 5/3, as indicated by the dotted lines. Right: Scaling of the 
relative density fluctuation Spjp, measured at the disc midplane. Saturation 
is reached for ~ 20 mG; below this value, the slope is approximately linear. 



field end the wavelength approaches the disc thickness (even at the 
base of the active layer). These geometric conditions apply to the 
case of ideal MHD, and hence are in addition to the requirement 
A > 1, describing the resistive quenching of the MRI. Note that be- 
cause of the z-dependence of the gas density p(z), the vertical band 
supporting MRI-unstable wave numbers will successively shift to- 
wards the midplane when increasing B-. In the limit of strong fields, 
this means that the MRI will be increasingly affected by the A > 1 
criterion. At the same time, the vertical extend of the MRI-active 
band situated between the magnetically-dominated halo and the 
resistively-quenched dead zone is expected to shrink, eventually 
shutting off the MRI altogether. 

Within the permissible range described above, we find that the 
Maxwell stress follows the scaling oc /l^^, which is considerably 
steeper than the linear dependence proposed by| Pessah et al.|p"007^ 
- see their fig. 2 - and somewhat steeper than the 3/2 dependence 
suggested by jSano et al.| ( |2004| l. We stress that this does not pose 



any inconsistency as the earlier results were obtained for the sim- 
plified case of unstratified MRI. |Pessah et al.]j2007^ , in fact, point 
out that their results are likely to be modified when accounting for 
the buoyant loss of magnetic energy in the stratified case (this loss 
is clearly seen in Figs.p^andfTTjbelow). 

In passing we note that the canonical mass accretion rate ob- 
served in T Tauri stars is 10"* Mq yr"' , and for disc masses similar 
to the minimum mass model this translates into a viscous stress pa- 
rameter ffss ~ O.OI ( [Hartmann et al.|1998| l. This value can indeed be 
achieved in models below the strong-field limit of the MRI. Setting 
aside the significant uncertainties in constraining T Tauri accretion 
rates, this raises the question whether or not a net field of around 
20 mG (required to generate the observed typical mass accretion 
rate) can be sustained by a self-consistent accretion-disc dynamo, 
or can be dragged in by the accretion flow that is fed by infall from 
the parent molecular cloud. 

While the turbulent stress, at very best, shows a weak break 
around field strengths of about 10 mG, the derived midplane den- 
sity fluctuation SpIp (see right-hand panel of Fig.|9| shows a clear 
saturation in the strong-field limit. To understand this saturation, 
we again need to take a look at the vertical structure of the disc. 



4.3 Disc structure at varying net-flux 

Let us recall that in Section |3] the density fluctuations at the disc 
midplane appeared to be determined by the wave energy available 
at the dead zone interface defined by the A = 1 condition. For mod- 
els Dl.l -4, the region where the MRI operated was limited by the 
column density - i.e. by the disc material being more efficient in 
shielding exterior sources of ionisation at higher mass loading. Be- 
cause the Elsasser number at the same time depends on the Alfven 
speed, and hence the magnetic field strength, we expect the position 
of the layer interface to depend on the external net flux. 

This dependence is clearly seen in Fig. [TO] where we plot the 
same quantities as in Fig. |2] along with the A = I interface for 
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Figure 10. Space-time evolution of the mean toroidal field (upper row), total Maxwell stress (middle row), and relative density fluctuation (bottom row) for 
run Dl-NVFa. The solid black line shows the transition between the MRI-active region and the dead zone, approximated by the condition A = x'^^HShf) = 1. 
Segmentation indicates intervals with dift'erent net vertical field strengths of B- = 2.68, 5.37, 10.73, 21.47 mG, respectively. The third interval, from 90 to 130 
orbits corresponds to model Dl.l, with a net field of 10.73 mG. When obtaining time-averaged quantities, the first 10 orbits of each interval are disregarded. 



the nm Dl-NVFa. We remind the reader that in this simulation we 
increase the net-flux by a factor of two every forty orbits to exam- 
ine how the dead zone structure changes (along with the stochastic 
forces experienced by embedded planetesimals). The half- width of 
the respective dead zones are 2.22±0.03, 1.89±0.04, and 1.67±0.04 
for the first three intervals indicated in Fig. [lO] This implies that 
doubling the net flux, in this regime, decreases the extent of the 
dead zone by about a third of a scale height - a trend that is equiva- 
lent to the one seen in Sect.[3]when reducing the column density by 
a factor of two. We note that this direct correspondence, however, 
is only valid for weak fields. 

To study stronger net vertical fields, we have to overcome an 
obstacle which becomes obvious in the last segment of Fig. [TO] 
where we observe a strong build-up of toroidal net flux caused by 
radial field. Diffusing into the midplane region, this radial flux gen- 
erates strong azimuthal fields through the differential rotation. Ac- 
cording to Kim & Ostriker (2000), the linear growth rate for ax- 
isymmetric MRI modes goes asymptotically to zero at low f}^ if 
the field has a significant toroidal component; this likely explains 
the disappearance of the field reversals and density streaks (asso- 
ciated with channel modes) as grows. As discussed in detail in 
Sect. |B2| the winding-up of azimuthal field is artificially enhanced 
within the shearing box approximation and would be counter-acted 
by radial diffusion within a global setup. While we thwart this ar- 
tificial build-up via an extra dissipative term l |B4[ l in the induction 
equation, the adopted time-scale turned out to be insufficient for 
net fields exceeding 20 mG. As a remedy, we devised an additional 
simulation run 'Dl-NVFV adopting a more conservative estimate 
for Tdiff, with dissipation faster by a factor of four. For the new run 
we chose an initial net field of 10.7 mG, which enables us to mon- 
itor the effect of the increased field dissipation term by comparing 



with the equivalent segment of run Dl-NVFb. As can be verified 
by inspection of the respective lines in Tab. [3] the derived quanti- 
ties agree to within 10% between the corresponding segments (i.e., 
orbits 100-130 of NVFa, and orbits 20-50 of NVFb). Moreover, 
the good agreement a posteriori supports our inclusion of the addi- 
tional magnetic diffusion term. 

Space-time plots for run Dl-NVFb are shown in Fig. [TT] 
where we see that for fields stronger than about 40 mG, MRI can- 
not be sustained and the disc returns to a non-turbulent state. In- 
terestingly, in this limit we still observe sporadic bursts of activity 
(partly restricted to one "hemisphere" of the disc). As above for run 
Dl-NVFa, we infer the half-thickness of the dead zone by means 
of the A = 1 criterion, and we yield 1.66+0.03, 1.54±0.04, and 
1.49±0.10//, for a net-flux field B- =10.7, 21.5, and 43 mG, re- 
spectively. We conclude that by increasing the field strength we 
narrow the MRI-active region sandwiched between the dead zone 
and magnetically-dominated halo, at the same time approaching a 
limit in the vertical position of the dead zone interface. Accord- 
ingly, the magnitudes of fluctuations Sp at the midplane saturate 
along with the stochastic torques, as displayed in the right and left 
panels of Figs.|9]and[T2] respectively. 

4.4 Gravitational torques as function of net-flux 

Unlike the runs with varying disc mass discussed earlier in this 
paper, we did not evolve the trajectories of embedded particles 
over time periods of ~ 200 orbits for each value of the magnetic 
field. Such an approach would be prohibitively expensive com- 
putationally. Instead, for most models, we monitor the stochastic 
torques experienced by embedded particles by accumulating time 
series, allowing us to infer orbital evolution through analysis of the 
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Right: Scaling with the density fluctuation Sp for runs Dl-WF (left-most 
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widths of the torque distributions o"(r,,), and the correlation times 
Tc (cf. Sect. |3.5| above). For model Dl-WF, corresponding to a weak 
field of 2.7 mG, we evolved the particle trajectories for 200 or- 
bits, and found excellent agreement between the simulation results 
and expectations based on the simple diffusion model discussed in 
Sect. |3.5[ thus justifying our general approach. The obtained values 
can be found in the eighth and ninth column of Tab. [3] respectively. 

As expected from the density fluctuations plotted in Fig.|9] the 
gravitational forcing, for magnetic field values below B, =^ 10 mG, 
depends strongly on the level of the applied field (see left-hand 
panel of Fig. |12[l. For field strengths exceeding 10 mG, however. 



we find saturation of the fluctuating torques, and this is consistent 
with the trend foimd in the density fluctuations in the preceding sec- 
tion. Moderately stronger fields are already MRI-stable and models 
in this regime are likely to lead to a chaotic time-behaviour of the 
disc as a whole. Given that we already reached accretion stresses 
compatible with limits given by observations of T Tauri discs for 
B; = 21.5 mG, we reckon that the occurrence of even stronger 
fields seems unlikely. In this sense, the obtained torques of around 
0.05 X 10** cm^ s"^ appear to be a firm upper limit for the given dead 
zone structure. Higher torques are hence only expected for higher 
ionising fluxes (cf. model D2 from Paper II) or lower dust abun- 
dances (not considered here, but see [Turner & Drake|2009^ , both 
of which decrease the height of the dead zone. 



To conclude our analysis of the results, we note that the range 
of values in 6p and cr(r,,) observed in the three runs Dl-WF, and 
Dl-NVFa/b can be used to infer the correlation between the two 
quantities. Recall that for models Dl.1-4, we found a somewhat 
surprising relation, where a higher value in absolute density fluc- 
tuations Sp would not translate into higher torques. In the limit of 
high net flux (and accordingly, high density fluctuations), this be- 
haviour is compatible with the scatter seen in the right-hand panel 
of Fig. [12] where we plot data points sampled from runs Dl-WF, 
and DI-NVFa/b. For weaker turbulence (i.e. for lower values of 
B,), however, we approximately recover the naturally expected lin- 
ear scaling of cr(r,.) with the density fluctuations. We conclude that 
only looking at the rms-value of Sp in general provides insufficient 
characterisation of the underlying gravitational stirring - a more 
sophisticated analysis is required to make quantitative predictions. 
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Figure 13. Absolute value of the rms density fluctuation Sp as a function 
of the vertical position for a model based on D1.4, but with Bq = 5.4 mG. 
Within the dead zone, excellent agreement with the theoretically predicted 
scaling based on the exponential background density profile is seen. 



5 THE RELATION BETWEEN DENSITY WAVE 
AMPLITUDE AND GRAVITATIONAL TORQUES 

One motivation of this work was to derive intuitive scaling rela- 
tions that would allow stochastic torque amplitudes to be expressed 
as simple functions of the disc model parameters. Such functions 
would provide straightforward prescriptions for gravitational stir- 
ring by turbulent discs with application to secular evolution models 
of planetesimal accretion, et cetera. The original strategy to achieve 
this aim was to break the problem down into three parts, i.e. (a) pre- 
dicting the extent of the dead zone and the amplitude of the turbu- 
lence as a function of disc mass and external net flux, (b) predicting 
the resulting amplitude of the induced spiral density waves (char- 
acterised by the rms density fluctuation Sp near the disc midplane), 
and (c) predicting the associated gravitational torque experienced 
by the embedded solids. The third step is particularly important as 
the evaluation of gravitational torques within direct simulations is 
computationally very expensive when large numbers of particles 
are being evolved. 

In a recent paper, [Okuzumi & Hirose| ( [2011^ have derived sim- 
ple predictor functions covering the first two steps outlined above. 
With respect to the third step, it is natural to assume that the result- 
ing torques scale linearly with the density fluctuation dp, but con- 
trary to our own expectations, we find that this assumption is not 
justified in the presence of a dead zone. Knowledge of dp alone is 
insufficient to construct a simple scaling relation that describes the 
gravitational stirring of planetesimals. Indeed the results presented 
in Sect.|3] are in stark conflict with this conjecture. In the following, 
we will elucidate how the discrepant findings can be reconciled. 
In our discussion, we are guided by the simple picture of density 
waves being excited in the active layer and then propagating into 
the dead zone, where they induce the gravitational stirring expe- 
rienced by the planetesimals. We will further refine this picture to 
include the effect of the background difl'erential rotation, as appears 
necessary to explain the results of our simulations. Finally, we de- 
fine a set of analytic formulae that provide a means of estimating 
stochastic torque amplitudes as a function of disc parameters. 

5.1 Effects due to background shear 

We recall from the discussion in Sect. |3.1| that absolute density fluc- 
tuations at the disc midplane should obey a scaling ^p^id sJPmid, 
which is indeed observed when comparing models Dl.l, D1.2, and 
D1.4, respectively. The agreement in predicting Sp is independently 
confirmed when looking at vertical profiles of Sp(z)- While we only 
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Figure 14. Model overview of the Sp versus torque relation. The shaded 
data points with error bars are taken from Fig. |12| and complemented by 
time-averages taken over ; = 30 orbits (diamonds). These points represent 
model Dl.l at varying external flux. Additional points (squares) are as la- 
belled. Dashed lines indicate a linear and constant scaling, respectively. 



used two particular positions (i.e. the interface and the midplane) 
in our previous discussion, we note that the derived relation should 
equally hold for any position within the disc, i.e., Sp(z) -y'p(z). 
This is plotted in Fig. [T3] where error bars indicate the deviation 
arising from temporal fluctuations within the adopted time interval 
? = 20 - 85 orbits. The excellent agreement with the predicted scal- 
ing confirms the assumption based on conservation of wave energy, 
and leads to the conclusion that step (b) in the procedure outline 
above is justified. 

We now contrast the assumption of linear dependence between 
the absolute density fluctuation Sp and its associated torque cr(r,.) 

- to be used in step (c) - with the results obtained from a compre- 
hensive set of numerical simulations. Figure [14] compiles all avail- 
able data points in the Sp versus 0"(r^) plane, including the ones 
previously shown in the right-hand panel of Fig. [T2] (i.e. for the 
dependence on the vertical net flux). We have also included addi- 
tional runs in Fig. [14] performed to provide improved coverage of 
the cr(r,.) - Sp plane. These runs were for the heavier disc D1.4 but 
with external magnetic field values B. = 5.4 and 2.7 mG (labelled 
in the figure). Clearly the scatter seen in the plot cannot be rec- 
onciled with a single linear relation between Sp and the resulting 
torque. Let us hence focus on a few particular sub-sets of models 
first. 

Evidently, the torque remains constanj^ when going from 
model DLl to DL2, and D 1.4 as discussed already in Section]?] In 
contrast, the torques for model Dl-NVFa/b (with varying net-flux 
at constant disc mass), are consistent overall with a linear scaling 

- although the torques exhibit a somewhat steeper than linear fall- 
olf as a function of Sp for Sp^ 0.1. The same trend holds for the 
additional set of runs based on model D1.4 but with varying net 
flux (cf . the red squares on the lower right in Fig. |14^ - here very 
little deviation from the linear scaling is seen. The most apparent 
discrepancy arises between model Dl.l and a four- times heavier 

* Note that a similar pattern (of a too-low torque at increased disc mass) 
holds when comparing Dl.l and D1.4 at weaker fields (i.e., 5.4mG, and 
2.7 mG), albeit the effect is seen strongest at 10.7 mG. 
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Figure 15. Top.' Contours of the two-point density correlation function for 
model D 1.2 at various vertical positions in the disc. The five dark contours 
show where the functions fall below 1/e of their peak values; light con- 
tours represent the first zero-crossings. Bottom: Comparison of the two- 
point ACF computed within z = ±H/2 of the disc midplane, for models 
Dl.1-4 (also cf. figs. 3 & 6 in Paper I), showing the same trend as above, 
but now due to the different width of the respective dead zones. 



disc at half the vertical net-flux, labelled 'D1.4 (5.4 mG)'. While 
these models agree in the midplane value of dp to within 10%, 
they differ in their effective torques by a factor of two. This clearly 
demonstrates that there has to exist a significant qualitative differ- 
ence in the density structure between the two models. The obvious 
difference between the strong-field/low-mass model and the weak- 
field/high-mass model is the width of the resulting dead zone. As 
we will see in the following discussion, this observation provides 
the key to understanding the apparent discrepancy. 



5.1.1 Two-point density correlation 

In the following, we build on the picture of linear sound waves 
transporting density perturbations into the dead zone region as out- 
lined in Sect. |3.1| We conjecture that the shearing-out of these ra- 
dially and vertically propagating spiral density waves affects the 
torque distribution inferred near the disc midplane. This picture 
gains support from looking at the two-point correlation of the gas 
density, which we compute in horizontal slabs of vertical extent 
z = ±H/2, and which we average over t = 20 - 220 orbits. The 
two-point ACF of the density field serves as a qualitative indicator 
for the degree of coherence in the density structures and was used in 
Paper I to compare local and global models of fully active MRI. To 
trace the characteristic evolution of the propagating density waves, 
we compute two-point ACFs at different vertical positions in the 
disc, starting from the dead/active zone interface and moving to- 
wards the midplane in steps of H/2. This is exemplified for model 
D 1 .2 in the upper panel of Fig. [15] but the observed trend holds for 



all the studied models. The plot clearly shows that density features, 
on average, are more sheared-out towards the midplane. 

We further conjecture that the stretching arises because the 
density waves become increasingly sheared out by the background 
differential rotation as they propagate from the vertical position 
where they are created down toward the midplane. To estimate the 
amount of stretching, we assume a vertical propagation speed of 
Cs = HQ., resulting in a propagation time Af = 2H/cs = 211 ' for 
a dead zone half-width of two scale heights. Assuming Keplerian 
shear with q = 3/2, and evaluating at a position x =^ 0.5 H away 
from the coordinate origin, this amounts to an azimuthal displace- 
ment of Ay = qQxAt = 3/4 HQ. X 2Q'' = 1.5 H, which is in 
decent agreement with the level of stretching seen in Fig. [15] As a 
further step, we checked that the characteristic pitch angle for den- 
sity waves observed at the dead zone interface is the same for all 
the models. 

With this knowledge at hand, we can now attempt to reconcile 
the results from Sect.jS] Given the significant widening of the dead 
zone when increasing the disc mass, we expect more sheared-out 
ACFs for model D1.4 than for model D1.2, and Dl.l, respectively. 
This is because the travel time for waves to reach the planetesimals 
situated near the disc midplane directly depends on the width of the 
dead zone. The suggested trend is clearly seen in the lower panel of 
Fig.[T5] where we plot the midplane ACFs for the three models. 

The shapes of the ACFs are characteristic of trailing spiral 
waves, and the degree of correlation in the azimuthal direction 
shows a clear trend with disc mass. This, in turn, supports the no- 
tion that the stronger fluctuations Sp in the heavier disc models 
are compensated by the more elongated aspect ratio in the density 
structures. It is well known that the winding-up of the spiral struc- 
tures due to differential rotation leads to a linear increase in the ra- 
dial wave-number /c, oc (3Q/2)tky in time. This process goes along 
with a phase-folding of the azimuthal structures. The enhanced az- 
imuthal symmetry of increasingly sheared out density waves, as 
viewed by embedded planetesimals, leads to a reduction in the av- 
erage amplitude of the stochastic torque due to partial cancellation 
of the induced azimuthal gravitational acceleration. 

5.1.2 Modified scaling relation 

So far, our reasoning has been largely qualitative. In this section we 
aim to provide more quantitative evidence of how the simple lin- 
ear relation between density fluctuations and gravitational torques 
might have to be modified. Having identified the distance away 
from the dead zone interface as the relevant parameter determining 
the pitch angle of the density field (and hence the effective gravita- 
tional force), it appears natural to use the dead zone half-width h^z 
as a parameter to empirically correct our prediction for the result- 
ing torques. This is done in Fig.[T6] where we boost the measured 
torques by a factor (/jdz/2//)'', with p to be determined. The rescal- 
ing is done to restore, as accurately as possible, a linear dependence 
a{Xy) «i 5p. The particular choice of the fiducial width 2H implies 
that we pivot around model D1.2, which has this dead zone width 
and hence remains unaffected. We point out that the rescaled torque 
amplitude does not have any physical implication but merely serves 
as a benchmark. In practice, one can use the inverse scaling factor 
to correct torque estimates based on predicted values of Sp for the 
attenuation caused by the stretching. By varying p and inferring 
deviations from a power-law, we find a best-fit index of p ^ 5/4, 
which does a reasonable job in collapsing the scatter plot to a sin- 
gle line in log-log space. As indicated by stars in Fig. [16] we ex- 
cluded the data points with field strengths higher than the nominal 
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Figure 16. Same as Figure [T4| but with the torques rescaled by (/idz/2//)'' 
to compensate for the attenuation of the measured torques by the stretching- 
out of the density field. Excluding the strong-field models (see text), we 
obtain a best-fit correction parameter of p ^ 5/4. 



value of 10.7 mG. This has been done because in this limit the dead 
zone width does not change any more, and the stirring amplitude 
saturates accordingly. We note, moreover, that the resulting slope 
is somewhat smaller than the expected linear relation. Given the 
simplicity of the approach, however, we regard this a satisfactory 
result. As can be seen in Fig. [16] the rescaling mostly affects mod- 
els with different column densities, and the models Dl.l, D1.2, and 
D1.4 now show the expected linear scaling with Sp. This also af- 
fects the heavy disc models with weaker fields (D1.4 [2.7 mG], and 
D1.4 [5.4mG]), which have a wider dead zone. We remark that the 
rescaling does not destroy the good fit already seen in these models, 
but at the expense of a power-law index deviating from one. More- 
over the slight deviation from the linear relation within the data 
points belonging to difl"erent time intervals in models Dl-NVFa/b 
(diamond-shaped points & Dl.l) is now much smaller. Overall we 
regard the procedure as reasonably successful in correcting for the 
attenuation of the torque amplitudes sufi'ered because the density 
waves shear out as they propagate towards the disc midplane. We 
hence suggest to compensate torque estimates based on (5pniid by 
an extra factor (hoz/2Hy'^^ with respect to a fiducial model with 
haz = 2H. Because heavier discs have wider dead zones and hence 
weaker effective torques (compared to the original scaling based 
on 5pmid alone), in practice, this leads to more favourable condi- 
tions for planetesimal accretion in heavier discs, where gas-drag 
damping of the induced velocity dispersion is enhanced. 



5.2 A simple stochastic torque prescription 

Having derived the additional correction factor for the stochastic 
torque amplitude as a function of hoz, we are now in the position 
to provide a simple analytic fitting formula that encapsulates all the 
data points of the studied set of models. We here focus on provid- 
ing rms torque amplitudes Ty and correlation times t^. Owing to 
saturation effects, we limit the scope of applicability to magnetic 
field strengths of B- ^ 10.7 mG. We also point out that the derived 
scalings are only strictly valid for the chosen ionisation model, rep- 
resentative of a a protosolar nebula at a distance of a = 5 au, and 
for a range of column densities around a few times the MMSN. 



As we have seen in sections[3]and|4] the width of the dead zone 
changed by about a third of a scale-height when either doubling 
the disc mass, or halving the magnetic net-flux, respectively. More 
precisely, we obtain a half-width 



hoz = 1.66/f 



A/lB 



A/is 



ln(S-/10.7mG) 
in2 

ln(2:/135gcm-^) 
in2 



(10) 



relative to model Dl, and with coefficients A/ij = 0.33 H, and 



Aha = 0.28// -0.04// 



ln(£/135gcm-^) 
ln2 ■ 



(11) 



i.e., accounting for a slightly weaker dependence of hoz on the net 
magnetic field when going to heavier disc models. 

The scaling of the rms midplane density fluctuation 6p with 
disc mass was derived in Sect. |3.1| Based on the empirical finding 
of i5pniid/pmid p^,y*^, we here use a power-law index of 0.37, i.e., 
deviating slightly from the analytically predicted square-root de- 
pendence. For the scaling with the vertical net flux, the assumption 
of a simple linear dependence between Sp and S; produced equally 
good results for the models based on Dl.l, and D1.4, respectively. 
Overall we thus arrive at a predictor function for the rms stochastic 
torque amplitude of 



r.v - r.v.Di 



135 gem 



B- 



10.7 mG / \ 1.66// 



/in 



(12) 



relative to the fiducial model Dl with r,,Di = 0.054 x ID** cm^ s"-, 
and with /?dz as given by \10\ , and (jTTJ. This formula ~ together 
with a canonical value of - 0.3 orbits - completes our simple 
prescription for including stochastic torques in generic planetesi- 
mal evolution models. 



6 LONG-TERM EVOLUTION 

In Sects. [33] and |4.4| we discussed the rate at which the velocity dis- 
persion of planetesimals is excited by gravitational interaction with 
turbulent density fluctuations in the disc midplane. The driving of 
the velocity dispersion is counteracted by the combined effects of 
gas drag and collisional damping pda et al.|2008^ , leading to a well 
defined equilibrium value for the rms random velocity. We now 
calculate the approximate magnitude of this equilibrium velocity 
dispersion and compare it with estimates for the catastrophic dis- 
ruption thresholds of planetesimals composed of strong or weak 
materials. 



6.1 Planetesimal equilibrium velocity dispersions 

Following closely the analysis presented in sect. 5.1 of Paper II, 
we estimate the equilibrium velocity dispersion, v^isp, for planetes- 
imals of different sizes, and compare it with the catastrophic dis- 
ruption thresholds for strong and weak aggregates to determine 
under which conditions collisions between (equal-sized) planetesi- 
mals will lead to growth rather than destruction. 

Balancing the rate of eccentricity excitation obtained in sim- 
ulations with the gas drag damping rate, we obtain the following 
expression for the resulting velocity dispersion (see Paper II for 
details) 



Vdi.sn 



IO'CdP 



1/3 



(13) 



© 2002 RAS, MNRAS 000,[T|j20l 



16 Gressel, Nelson & Turner 



Rubble piles 



Collislonol 
disruption 
regime 



- Strong 

- oggregotes 
survival 



: Weok 

- oggregotes 

- survivol 
regime 




• No 


deod zone 


• 1 X 


MMSN 


• 2 X 


MMSN 


-•- 4 X 


MMSN 


4 X 


MMSN (5.4mG) 




10.00 


100.00 



1.00 

plonetesimal size R [km] 



_ Collislonol 
disruption 
regime 



Rubble piles 




Weok 

oggregotes 

survivol 

regime 



* B = 


43.0mG 


* B = 


21.5mG 


• B = 


10.7mG 


* B = 


5.37mG 


• B = 


2.68mG (scoled) 


B = 


2.68mG (direct sim.) 




10.00 


100.00 



I 



1 .00 

plonetesimal size R [knn] 



Figure 17. Left: Variation of the equilibrium velocity dispersion as a function of planetesimal size and the disc mass. Increasing the disc mass increases 
the midplane density and planetesimal collision probability, thereby increasing gas drag and coUisional damping. Right: Variation of equilibrium velocity 
dispersion as function of planetesimal size and magnetic field strength. 



Similarly, balancing the rate of eccentricity growth with the 
mutual collisional damping rate for planetesimals (assuming all 
solids in the disc have grown into planetesimals of a given size, R^, 
with a mean density Pp = 2 g cm"', and the coefficient of restitution 
- 0) gives the expression 



Vdis 



(14) 



The catastrophic disruption thresholds for strong (basalt) and weak 
aggregates are given by Stewart & Leinhardt| |2009 ), and are dis- 
cussed in detail in Paper II in the context of planetesimal stirring 
in turbulent discs. We plot the catastrophic collision velocities in 
Fig.[T7]as a function of planetesimal radius (solid black lines). We 
assume that the largest 100 km sized bodies are rubble-piles rather 
than monolithic bodies, and so adopt the appropriate disruption ve- 
locity in this limit. In this figure we also plot the equilibrium ve- 
locity dispersions for planetesimals of different size from the dif- 
ferent simulations. The plotted values correspond to the smaller of 
those calculated using eqns. (jTSj and l |14^ . The left panel shows re- 
sults for different disc masses with a fixed external magnetic field 
strength of 10.7 mG. The right panel compares results for different 
magnetic field strengths in a disc with fixed mass (approximately 
equal to the MMSN model). 

The left panel of Fig. ^] shows that planetesimals with radii 
between 0. 1 and 1 km have equilibrium velocity dispersions that 
lie between the catastrophic disruption thresholds for weak and 
strong materials when embedded in discs with masses between 1 
and 4 times the MMSN. Collision outcomes for these bodies will 
depend sensitively on the material properties. 10 km bodies in discs 
with 2 or 4 times the MMSN mass will have velocity dispersions 
that are marginally lower than the catastrophic disruption thresh- 
olds for even weak materials. lOOkm-sized rubble piles are clearly 
safe from disruption in all models. The decreasing equilibrium ve- 
locity with increasing disc mass arises because the gas drag scales 
with the midplane density, the mutual collision rate depends on the 
number density of planetesimals, and the strength of the stochastic 
stirring was found to be essentially independent of the disc mass 
for these models. 



Examining the right panel of Fig. [TT] we see that reducing 
the value of the magnetic field leads to a reduction in the stiiTing 
of planetesimals. Decreasing the field strength increases the size 
of the dead zone and reduces the level of turbulence in the active 
layers. We remind the reader that for most simulations conducted 
with varying magnetic field strength, the orbital trajectories of par- 
ticles were not evolved under the influence of disc forces, but in- 
stead we monitored the time varying torque due to the turbulent 
fluctuations acting on the particle swarm. Estimates for the values 
of Ccr(e) used in eqns. ^13\ and ^14\ to calculate the equilibrium 
velocity dispersions were obtained by assuming a linear scaling 
between the rms fluctuating torque and C^(e), using the run with 
B- = lOmG to provide the normalisation. The validity of this ap- 
proach is confirmed by the run with B- = 2.68 mG for which parti- 
cle trajectories were evolved under the influence of disc forces, and 
for which both the predicted and simulated values of Ca-{e) were 
found to be in very good agreement, as shown in the right panel 
of Figure [T7] For this run, which contains the weakest field con- 
sidered, the equilibrium values of I'jjjp fall beneath the catastrophic 
values appropriate for weak aggregates for all planetesimal sizes, 
suggesting that all bodies will be safe from collisional destruction 
in such a disc. For stronger fields, however, the values typically lie 
between the thresholds for strong and weak materials. lOOkm-size 
rubble-piles are again safe from collisional break-up in all models. 
We note that the predicted mass accretion rate for the disc model 
with = 2.68 mG is rather low, being M ^ 5 x 10"' Mq yr"' . This 
is a factor of two lower than the canonical value of 10"* Mq yr"' 
observed for T Tauri stars (Sicilia-Aguilar et al.|2004^ . The scaling 
we obtain with field strength and disc mass, however, suggest that a 
slightly stronger field can generate a larger accretion rate in the ac- 
tive zone, and in a more massive disc than the MMSN the velocity 
dispersion will remain lower than the catastrophic values for weak 
materials. This expectation is supported by the simulation D1.4b 
described in Tab. [3] which used the heaviest disc model DI.4 and 
an external magnetic field S; = 5.4 mG. Turbulent stresses induce a 
mass accretion rate ~ 2.7 x 10"* Mq yr"' , and the values of I'^isp are 
illustrated by the lowest set of points in the left panel of Fig.[T7] It is 
clear that disc models can be constructed that satisfy requirements 
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on disc accretion rates and allow even the weakest planetesimals to 
avoid collisional disruption. It is advisable not to over-interpret the 
results of our study, however, due to the simplifying assumptions 
that have gone into our calculations. Suffice it to say that we have 
been able to produce disc models with reasonable accretion rates 
within which the dead zone provides a safe-haven against catas- 
trophic disruption of planetesimals of all sizes even when they are 
composed of essentially strength-less materials. 

6.2 Implications for planet formation 

We now summarise the key results from our studies of planetesi- 
mals embedded in turbulent discs, focusing on their relevance for 
planet formation theory. We do this by first providing a brief review 
of our findings from Papers I and II and then discuss how these are 
affected in view of the new results obtained in this paper. 

In Paper I, we considered the dynamical evolution of small 
particles (pebbles and boulders) whose interaction with the gas was 
dominated by gas drag, in addition to larger bodies (planetesimals) 
whose interaction is dominated by gravitational interaction with 
turbulent density fluctuations. Particles tightly coupled via gas drag 
have sizes < 1 m at 5 au, and were shown to rapidly achieve a veloc- 
ity dispersion equal to the gas turbulent velocity. At the midplane in 
a dead zone this velocity is typically 10-30 m s"' , possibly slowing 
protoplanet growth in models that rely on accretion of small parti- 
cles by planetesimals and embryos (e.g. |Rafikov|2004||Johansen &| 
|Lacerda|2010l|Ormel & Klahr|2010| l. 

We now focus on the evolution of larger bodies whose dynam- 
ical evolution is determined by gravitational interaction with turbu- 
lent density features. Papers I and II demonstrate that a model of 
planetesimal formation based on gradual binary sticking of smaller 
particles cannot operate in a fully turbulent disc due to rapid growth 
of the velocity dispersion above catastrophic disruption thresholds 
for bodies of size < 10 km. Models that invoke rapid formation 
of large (i.e., ~ 100km) planetesimals may operate (e.g. |Johansen| 
let al.|2007[[Cuzzi et al.|2"008} , but runaway growth to form larger 
bodies is not possible because the velocity dispersion quickly rises 
to the surface-escape velocity of these objects. Forming planets 
within the disc lifetime appears to be difficult in a fully turbulent 
disc without a dead zone, unless a model involving gravitational 
instability is invoked ( |Boss|1998] l. 

Results from Paper II (and model Dl.l in this work) show 
that stochastic forcing of planetesimals is ~ 20 times weaker for a 
minimum-mass disc with a nominal dead zone than in a disc with- 
out a dead zone. Equilibrium velocity dispersions for 100 m - 10 km 
bodies lie between disruption thresholds for weak and strong ma- 
terials ^Stewart & Leinhar dt 2009), being between a few metres 
and a few tens of metres per second (cf. sect. 5.3 in Paper II for 
details). In principle, planetesimals that form via gradual accumu- 
lation can survive mutual collisions in the nominal model if com- 
posed of moderate-strength materials. To form in the first place, 
however, they need to overcome barriers to growth occurring at mm 
and metre sizes ( [Brauer et al.|2008l|Zsom et al.|20I0) . lOkm-sized 
bodies forming gradually will not experience runaway growth to 
form larger oligarchs as their velocity dispersion reaches ~ 10 m s"' 
within ~ 10^ yr. Because the dispersion remains below the destruc- 
tion threshold, further growth can proceed; eventually delayed run- 
away growth may occur for resulting bodies with ~ 100 km radii 
(also see |Ormel et al.|2010| l. 

Rapid formation of large planetesimals with sizes > 100 km 
through turbulent concentration of chondrules jCuzzi et al. '2008) or 
the streaming instability (Johansen et al.|2007] l in the nominal dead 



zone results in runaway growth because of slow turbulent stiiTing 
of the velocity dispersion. Self-stirring by planetesimals determines 
the dispersion in this case rather than turbulence. We note that all 
dead zone models predict that 100 km rubble-piles are safe against 
turbulence-induced collisional destruction. 

Increasing the disc mass by a factor of two or four does not 
dramatically change the above picture, as demonstrated by models 
D1.2 and D1.4 in Sect. [3] These models result in smaller equilib- 
rium dispersion velocities for larger disc masses because of the con- 
stant stirring strength and increased eccentricity damping. Increas- 
ing the external magnetic field strength above 10.7 mG increases 
stochastic forcing and vjisp, but saturation is attained when the ex- 
ternal field approaches values too strong for the MRI to operate 
in the active zone. Reducing the field strength weakens stochastic 
forcing by increasing the dead zone width and reducing activity 
levels in the turbulent layers. A model with B- = 2.7 mG results 
in Vdisp values below the catastrophic disruption threshold for the 
weakest aggregates. This important result demonstrates that weak 
100 m - 100 km planetesimals can survive collisional destruction in 
a reasonable dead zone model, albeit one whose mass accretion rate 
is a factor of two below the canonical mass accretion rate in T Tauri 
stars (but still within the observed range of values). Model D1.4b 
demonstrates, however, that a heavier disc with a slightly stronger 
field can generate a larger mass accretion rate and similar velocity 
dispersion (see Tab.[3]and Fig.[T7j. Weaker stochastic forcing leads 
to slower growth of the velocity dispersion toward the escape ve- 
locity of lOkm-sized bodies. The time for the velocity dispersion 
to reach 10 ms"' in the weak magnetic field model is ~ 3Myr, so 
runaway growth of lOkm-sized bodies can occur in this disc. 



7 CONCLUSIONS 

This is the third paper in a series that examines the dynamics of 
planetesimals embedded in turbulent protoplanetary discs. In the 
present work we have used local disc models situated at 5 au that 
include a simple equilibrium gas-grain chemistry to examine how 
changing the disc column density and the external magnetic field 
strength modifies the size of the dead zone, and the level of gravita- 
tional stirring of particles located at the disc midplane. We consider 
column densities that lie approximately between one and four times 
the minimum mass solar nebula values, and external magnetic field 
strengths between 2.68 - 43 mG. Our main conclusions can be sum- 
marised as follows: 

(i) Increasing the disc mass by a factor of two or four above 
the nominal value for the minimum-mass model leads to the un- 
expected result that the strength of particle stirring is independent 
of disc mass. Scaling relations based on a simple picture of linear 
waves being excited in the active layer and propagating into the 
dead zone predict that absolute density fluctuations at the midplane 
will scale with the square-root of the column density. This scaling 
is observed in the simulations, but gravitational stirring is found to 
be insensitive to the larger absolute density fluctuations. 

(ii) Analysis of two-point correlation functions for the midplane 
density fields between discs with dilferent dead zone heights show 
that typical density structures in discs with larger dead zones are 
more elongated in the azimuthal direction. Density waves propa- 
gating from the active layer into the dead zone are more sheared 
out when they reach the midplane in larger dead zones, resulting in 
reduced stochastic gravitational forces. 

(iii) We suggest this phenomenon explains the insensitivity of 
planetesimal stirring to increasing disc mass, as is demonstrated by 
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a scaling-law for the stochastic torques that depends on the dead 
zone vertical extent, and which accounts reasonably for our simu- 
lation results. 

(iv) Decrementing the external magnetic field strength increases 
the size of the dead zone and decreases the level of activity in the 
turbulent layers. This, in turn, leads to reduced absolute density 
fluctuations in the disc midplane and corresponding reductions in 
the level of gravitational stirring experienced by embedded plan- 
etesimals. 

(v) Moderate increments in the external magnetic field strength 
reduce the size of the dead zone, and increase the level of activity 
in the turbulent layers. This, at first, results in an increased stirring 
of planetesimals, but quenching occurs when the position of the 
dead zone interface reaches saturation. This is because the amount 
to which the active layers can constrict the dead zone is limited (via 
the A = 1 criterion) by the steep rise in rjiz)- For our model, torque 
fluctuations saturate for B- i 20 mG, while turbulent stresses grow 
up to net fields of ^ 40 mG, where the MRI eventually shuts off. 

(vi) A disc hosting a modest external field (of = 2.7-5.4 mG) 
leads to equilibrium velocity dispersions for 100 m to 10 km plan- 
etesimals below the catastrophic disruption threshold for even the 
weakest aggregates. This demonstrates that reasonable disc models 
can be constructed in which the dead zone provides a safe-haven 
for weak planetesimals against collisional destruction. 

There remain a number of caveats that we have not yet ad- 
dressed, and which may modify the detailed conclusions drawn in 
this paper. These include adoption of a simple equation of state 
for the disc, and a mono-disperse size distribution for planetesi- 
mals when estimating collisional damping rates. The most impor- 
tant omissions are the Hall effect and ambipolar diffusion from the 
induction equation, which may affect the size of the dead zone 
and the level of turbulent activity in the active zones ( |Wardle &| 
|Salmeron|2011| ). Nonetheless, we expect that the main result of our 
study remains robust, namely that dead zones lead to much reduced 
stochastic forcing of planetesimal random velocities, and provide 
potentially benign locations for planetary formation. 
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APPENDIX A: NUMERICAL CONVERGENCE CHECK 

To verify that our simulations are numerically converged, we per- 
formed two additional runs identical to model Dl.l, but with the 
numerical resolution changed to 16, and 32 grid cells per H, re- 
spectively. Owing to constraints in computational cost, we were 
not able to evaluate gravitational torques for the higher-resolved 
run but have to rely on the hydrodynamic quantities as proxies in- 
stead. The dependence of key quantities defining the turbulent state 
of the fluid are shown in Figure [AT] 

As was briefly discussed in Sect. |4.2[ the numerical require- 



ments for resolving the development of the MRI in the presence of 
a net vertical field are safely established by linear theory. Accord- 
ingly, we see very little change in the turbulent stresses ffMaxw, and 
Q'Rcyn when increasing the numerical resolution beyond the basic re- 
quirement. However, from our study of the box-size dependence in 
Paper I, we recall that convergence in these quantities does not nec- 
essarily guarantee that secondary effects are equally well resolved. 
In particular, this applied to the convergence (with box size) in the 
amplitude dp of the excited of spiral density waves. Because of the 
predicted non-linear steepening of the analytic wave profile l |Heine-| 
|maim & Papaloizou|201 1) , we suspect that spurious numerical dis- 
sipation might have an influence on the lifetime of wave fronts and 
hence on the overall amplitude of density perturbations. 

This effect is seen in the dependence of 5p on the numeri- 
cal resolution as shown by the upper line in Figure |A1| There 6p 
changes by almost 70% when increasing the number of points from 
16 to 24 per H, i.e., for a 50% increase in linear resolution. In con- 
trast, when going from 24 to 32 cells (corresponding to an addi- 
tional 33% increase in resolution), the change in 6p remains below 
15%, such that the discrepancy between the two resolutions is con- 
sistent within the error bars. Put into proportion with the level of 
accuracy aimed for, this indicates that our standard resolution of 
24/// is reasonably converged and will allow to predict turbulent 
density fluctuations to within about ±25%. 



APPENDIX B: MODIFICATIONS TO THE INDUCTION 
EQUATION 

Let us first focus on the time step limitations imposed by the high 
value of the diffusivity coefficient t] near the disc midplane (cf. 
Fig.[T|(. Given the need to evolve the simulations over several hun- 
dred dynamical time scales to follow the evolution of the embed- 
ded planetesimals, using these profiles is prohibitive. By inspect- 
ing earlier simulations and test runs, however, we have determined 
that there exists no significant small-scale structure in the regions 
of high diffusivity. This is hardly surprising as the diffusion term 
is very efficient at removing features in the magnetic field topol- 
ogy. We have found that imposing a cap to the diffusivity, 77, corre- 
sponding to a value of Rm = 3, does not significantly change the 
observed field topology, while at the same time permitting a much 
larger computational time step. However, as we will see below, ad- 
ditional care has to be taken to make this approach viable. 

Bl Evolution of mean fields 

Owing to the periodic nature of the shearing box approximation, the 
only field structures that remain in the dead zone region are verti- 
cal variations in the radial and azimuthal field. Magnetic fields are 
essentially uniform in the horizontal directions. The vertical vari- 
ations are given by the horizontally-averaged mean values B^{z), 
and By{z)- Let us now examine how these are aff'ected by diffusion. 
Note that the induction equation 

a,B- Vx(vxB-;7VxB) = 0, 

VB = 0, (Bl) 

is linear in both B, and 77, and can be decomposed accordingly. For 
reasons that will become obvious in a moment, we split the induc- 
tion equation into two parts, reflecting the evolution of the horizon- 
tally averaged mean fields and their fluctuations (cf . |Gressel|20 1 Q) . 
Because the geometric averaging operation trivially commutes with 
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20 Gressel, Nelson & Turner 



the differential operators in ([BTJ, this decomposition is exact. Fo- 
cusing on the evolution of the mean fields, we now write down the 
diffusive part of the related induction equation: 

d,B,{z) = ■■■ + d,[ niz) d-Mz) ] , ie[x,y} . (B2) 

Introducing a cap on the value of rj (and hence removing its z- 
dependence beyond a certain depth within the dead zone) may sig- 
nificantly alter the character of |B2](, and we illustrate this by means 
of a simple example: Note that for constant rj, a field varying lin- 
early with z is not subject to any further diffusion as rid\Biiz) = in 
this case. In contrast, any "bump" in r] will make rfiz) d-Bi(z) vary 
in z even if d^Bj(z) is constant. The diffusion term will then act to 
flatten a given linear ramp, removing any residual gradients in B(z). 
Looking at the strong gradients d. r]{z) seen in Fig.[T| this implies a 
potentially efficient means to redistribute field. Ultimately, a weak 
gradient in B can lead to strong diffusion in the presence of a large 
gradient in rj. 

Because it is a one-dimensional equation, the overhead of 
solving \B2\ is computationally inexpensive, even under stringent 
time step constraints. Consequently, we split ||BTJ into four parts 
by separating both B = B + B' and rj = ;7>3(x, t) + ti^t,(z, t). Here 
subscripts refer to the resulting magnetic Reynolds number, Rm, 
as indicated by a dashed line in Fig. [T| also note the function ar- 
guments indicating that the highly diffusive low-Rm part of (i.e. 
77<3) is only applied to the mean-field equation. In summary, this 
approach allows us to retain the full ?/ profile in |B2](, while using 
the capped i] profile in the expensive three-dimensional equation. 
We finally remark that retaining the full 77 in the mean fields will 
likely enhance the leakage of field into the dead zone, as discussed 
above. 



even simpler approach and estimate the time scale on which diffu- 
sion redistributes azimuthal flux as 

Tdiff = . where we chose = njr . (B3) 
kjrj 

Noting that a Fourier mode of wave number would exponentially 
decay according to exp(-?/Tdifr), we remove toroidal (and radial) 
flux via an additional source term 

d,B,{z) = .■■-B,(z)lTm (B4) 

in the one-dimensional part of the induction equation. As we will 
find when discussing the results, this choice of the global decay 
time scale retains some weak fields within the midplane (also cf. 
fig. 3 in |Hirose & Tumer||201 1^ , indicating that we are not sig- 
nificantly over-damping the system. As can be seen in the upper 
panel of fig. 5 of Paper II, model Dl is affected by relatively strong 
fields diffusing into the midplane region. To make the comparison 
with models D1.2, and D 1.4 as straight-forward as possible, we de- 
cided to re-run this simulation including the diffusive sink term in 
Eqn. (|B4j, giving the new version the label 'Dl.l'. 

This paper has been typeset from a Tj5f/ KTgX file prepared by the 
author. 



B2 Accounting for super-box scale effects 

Neglect of the large-scale radial structure makes the shearing-box 
a poor set-up for global magnetic field evolution. This deficiency 
becomes most apparent within the dead zone, as we have discussed 
briefly above. Accounting for the global disc topology, [Turner &| 
|Sano| poo s'! derived a criterion for the evolution of large-scale 
magnetic fields within dead zones. Owing to the limitations im- 
posed by the neglect of curvature terms, shearing boxes cannot re- 
produce the related threshold, which applies in global cylindrical 
coordinates. This is because the divergence of radial field lines is 
not included, so 6, can be uniform where 6,. would decline as 1/r. 
Omitting radial gradients in B, contributes to the uniform shear- 
generated B^ and the weak resistive diffusion. 

Moreover, in a global disc model, the Keplerian rotation pro- 
file Q.ir) oc r^^l'^ leads to a radial variation of the shearing time 
scale. In contrast, a single orbital frequency Qq applies throughout 
the local box, such that a uniform radial field stretches out into a 
uniform toroidal field. With the toroidal field uniform, the resis- 
tive term 77 B^ in the induction equation toroidal component is 
zero, preventing the diffusion of the shear-generated fields. As a 
result, local simulations of dead zones can exhibit the growth of 
large toroidal flux (cf. |Hirose & Turner|201l| , which may spuri- 
ously feed-back into the MRI-active regions. 

One possibility to account for the global field evolution would 
be to remove toroidal flux at a rate consistent with the solution 
of the global induction equation. For example, if one assumed 
a power-law radial dependence for the dead zone toroidal field 
B^ oc r'', and a radially uniform resistivity, the axisymmetric cylin- 
drical induction equation has resistive term c),i?0 = T] B^{p^ -i) r"^ . 
Because of the uncertainty related to the parameter p, we chose an 
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